home *** CD-ROM | disk | FTP | other *** search
/ Eagles Nest BBS 8 / Eagles_Nest_Mac_Collection_Disc_8.TOAST / Developer Tools⁄Additions / MacScheme20 / Mathlib / carith.scm < prev    next >
Encoding:
Text File  |  1989-04-27  |  16.8 KB  |  634 lines  |  [TEXT/????]
  1. ;;; $Header: carith.scm,v 1.5 88/09/01 22:13:23 GMT cph Exp $
  2. ;;;; Complex Arithmetic with Operators
  3.  
  4. ;;;; This file is patched by Hal to get rid of operators, making them the
  5. ;;;; same as 1-argument procedures.  I have done this keeping in stubs, so
  6. ;;;; it will be easier to change when we figure out the right thing.  I have 
  7. ;;;; NOT made the corresponding changes to GARITH.
  8. ;;;; I have NOT worried about conditionalizing this so that it works with
  9. ;;;; other than the MIT system.
  10.  
  11. (if-mit
  12.  (declare (usual-integrations = + - * /
  13.                               zero? 1+ -1+
  14.                       ;; truncate round floor ceiling
  15.                   sqrt exp log sin cos)))
  16.  
  17. ;;; Requires macro file <mumble>SYN.SCM (in GENERAL)
  18.  
  19. ;;; The definitions follow the order in
  20. ;;; Revised^3 Report on the Algorithmic Language Scheme.
  21.  
  22. (if-mit
  23.  
  24. (define-integrable r-complex? (access complex? system-global-environment))
  25. (define-integrable r-number? (access number? system-global-environment))
  26.  
  27. (define-integrable r-zero? (primitive zero?))
  28.  
  29. (define-integrable r-= (primitive &=))
  30.  
  31. (define-integrable r-1+ (primitive 1+))
  32. (define-integrable r--1+ (primitive -1+))
  33.  
  34. (define-integrable r-+ (primitive &+))
  35. (define-integrable r-- (primitive &-))
  36. (define-integrable r-* (primitive &*))
  37. (define-integrable r-/ (primitive &/))
  38.  
  39. (define-integrable r-sin (primitive sin))
  40. (define-integrable r-cos (primitive cos))
  41.  
  42. (define-integrable r-sqrt (primitive sqrt))
  43. (define-integrable r-exp (primitive exp))
  44. (define-integrable r-log (primitive log))
  45.  
  46. (define r-acos (access acos system-global-environment))
  47. (define r-asin (access asin system-global-environment))
  48. (define r-atan (access atan system-global-environment))
  49. (define r-expt (access expt system-global-environment))
  50. (define r-abs (access abs system-global-environment))
  51.  
  52. ) ;; End of if-mit
  53.  
  54. ;;; Real primitives are inherited from underlying Scheme.
  55.  
  56. (if-pcs
  57.  
  58. (move-definition complex? r-complex?)
  59. (move-definition number? r-number?)
  60.  
  61. (move-definition zero? r-zero?)
  62.  
  63. (move-definition = r-=)
  64.  
  65. (move-definition 1+ r-1+)
  66. (move-definition -1+ r--1+)
  67.  
  68. (move-definition + r-+)
  69. (move-definition - r--)
  70. (move-definition * r-*)
  71. (move-definition / r-/)
  72.  
  73. (move-definition sin r-sin)
  74. (move-definition cos r-cos)
  75. (move-definition atan r-atan)
  76. (move-definition acos r-acos)
  77. (move-definition asin r-asin)
  78.  
  79. (move-definition sqrt r-sqrt)
  80. (move-definition exp r-exp)
  81. (move-definition log r-log)
  82. (move-definition expt r-expt)
  83.  
  84. (move-definition abs r-abs)
  85. ) ;; End of if-pcs
  86.  
  87. ;;; Operator stuff
  88. ;;; For this version, "operators" are 1-argument procedures
  89.  
  90. ;;(define-integrable operator-type '*operator*)
  91.  
  92. ;;;we keep these as stubs, just in case we want to try something else.
  93.  
  94. (define-integrable &procedure->operator
  95.   (lambda (procedure) procedure))
  96.  
  97.  
  98. (define-integrable &operator?
  99.   (lambda (object) (procedure? object)))
  100.  
  101. (define-integrable &operator->procedure
  102.   (lambda (object) object))
  103.  
  104. (define identity-operator
  105.   (&procedure->operator (lambda (object) object)))
  106.  
  107. (define (operation*operator->operator operation operator)
  108.   (&procedure->operator
  109.    (lambda (object)
  110.      (operation ((&operator->procedure operator) object)))))
  111.  
  112. (define (operation*constant*operator->operator operation constant operator)
  113.   (&procedure->operator
  114.    (lambda (object)
  115.      (operation constant ((&operator->procedure operator) object)))))
  116.  
  117. (define (operation*operator*constant->operator operation operator constant)
  118.   (&procedure->operator
  119.    (lambda (object)
  120.      (operation ((&operator->procedure operator) object) constant))))
  121.  
  122. (define (operation*operator*operator->operator operation operator1 operator2)
  123.   (&procedure->operator
  124.    (lambda (object)
  125.      (operation ((&operator->procedure operator1) object)
  126.         ((&operator->procedure operator2) object)))))
  127.  
  128. (define (procedure->operator procedure)
  129.   (if (procedure? procedure)
  130.       (&procedure->operator procedure)
  131.       (error "Illegal operator definition" procedure)))
  132.  
  133. (define (operator->procedure operator)
  134.   (if (&operator? operator)
  135.       (&operator->procedure operator)
  136.       (lambda (object) operator)))
  137.  
  138. (define (operator-compose operator1 operator2)
  139.   (&procedure->operator
  140.    (lambda (object)
  141.      ((operator->procedure operator1)
  142.       ((operator->procedure operator2)
  143.        object)))))
  144.  
  145.  
  146. (define (operator-value operator object)
  147.   (if (&operator? operator)
  148.       (with-reasonable-object object (&operator->procedure operator))
  149.       operator))
  150.  
  151. ;;;this is the advertized name, in a world without operators
  152.  
  153. (define (limiting-value procedure object)
  154.   (operator-value (procedure->operator procedure) object))
  155.  
  156. (define (with-reasonable-object object procedure)
  157.   (fluid-let ((operator-value-retry
  158.                 (call-with-current-continuation
  159.                   (lambda (k)
  160.                     (lambda ()
  161.                       (set! object (perturb-object object))
  162.                       (k k))))))
  163.     (procedure object)))
  164.  
  165. ;;; Because of the following definitions, these operators are
  166. ;;; restricted to numbers.
  167.  
  168. (define operator-perturbation-epsilon 1e-11)
  169.  
  170. (define (perturb-object object)
  171.   (+ object operator-perturbation-epsilon))
  172.  
  173.  
  174. (define (operation*operator*object->operator operation operator object)
  175.   (cond ((or (real? object) (&complex? object))
  176.      (operation*operator*constant->operator operation operator object))
  177.     ((&operator? object)
  178.      (operation*operator*operator->operator operation operator object))
  179.     (else
  180.      (error "Bad number" object))))
  181.  
  182. (define (operation*object*operator->operator operation object operator)
  183.   (cond ((or (real? object) (&complex? object))
  184.      (operation*constant*operator->operator operation object operator))
  185.     ((&operator? object)
  186.      (operation*operator*operator->operator operation object operator))
  187.     (else
  188.      (error "Bad number" object))))
  189.  
  190. (if-mit
  191.  
  192. (define-integrable complex-type 60)  ;(microcode-type 'COMPLEX) not integrated
  193.  
  194. (define-integrable (&complex? z)
  195.   (object-type? complex-type z))
  196.  
  197. (define-integrable (&complex re im)
  198.   (system-pair-cons complex-type re im))
  199.  
  200. (define-integrable &real-part system-pair-car)
  201.  
  202. (define-integrable &imag-part system-pair-cdr)
  203.  
  204. ) ;; End of if-mit
  205.  
  206. (if-pcs
  207.  
  208. (define-integrable complex-type '*rect*)
  209.  
  210. (define-integrable &complex?
  211.   (lambda (z)
  212.     (and (pair? z) (eq? (car z) complex-type))))
  213.  
  214. (define-integrable &complex
  215.   (lambda (re im)
  216.     (list* complex-type re im)))
  217.  
  218. (define-integrable &real-part cadr)
  219.  
  220. (define-integrable &imag-part cddr)
  221.  
  222. ) ;; End of if-pcs
  223.  
  224. ;;;; Complex number data abstraction
  225.  
  226. (define (&conjugate z)
  227.   (&complex (&real-part z)
  228.         (r-- 0 (&imag-part z))))
  229.  
  230. (define (&sqr-magnitude z)
  231.   (define r-sqr (lambda (x) (r-* x x)))
  232.   (r-+ (r-sqr (&real-part z))
  233.        (r-sqr (&imag-part z))))
  234.  
  235. (define (&magnitude z)
  236.   (r-sqrt (&sqr-magnitude z)))
  237.  
  238. (define (&angle z m)
  239.   (if (r-zero? m)
  240.       (if (fluid-bound? operator-value-retry)
  241.           ((fluid operator-value-retry))
  242.           (error "angle: magnitude is zero"))
  243.       (r-atan (&imag-part z) (&real-part z))))
  244.  
  245.  
  246. ;;; Complex stuff
  247.  
  248. (define (&make-rectangular re im)
  249.   (if (r-zero? im)
  250.       re
  251.       (&complex re im)))
  252.  
  253. (define (&make-polar r theta)
  254.   (if (r-zero? r)
  255.       0
  256.       (&make-rectangular (r-* r (r-cos theta))
  257.              (r-* r (r-sin theta)))))
  258.  
  259. (define i (&make-rectangular 0 1))
  260. (define -i (&make-rectangular 0 -1))
  261.  
  262. ;;;; Bottom layer
  263.  
  264. (define-integrable (make-binary-componentwise-operation r-opr r-opi
  265.                             operator-ok?
  266.                             accum error-string)
  267.   (letrec ((operation
  268.         (lambda (z1 z2)
  269.           (cond ((real? z1)
  270.              (cond ((real? z2) (r-opr z1 z2))
  271.                ((&complex? z2)
  272.                 (accum (r-opr z1 (&real-part z2))
  273.                    (r-opi 0 (&imag-part z2))))
  274.                ((and operator-ok? (&operator? z2))
  275.                 (operation*constant*operator->operator operation
  276.                                    z1 z2))
  277.                (else (error error-string z2))))
  278.             ((&complex? z1)
  279.              (cond ((real? z2)
  280.                 (accum (r-opr (&real-part z1) z2)
  281.                    (r-opi (&imag-part z1) 0)))
  282.                ((&complex? z2)
  283.                 (accum (r-opr (&real-part z1) (&real-part z2))
  284.                    (r-opi (&imag-part z1) (&imag-part z2))))
  285.                ((and operator-ok? (&operator? z2))
  286.                 (operation*constant*operator->operator operation
  287.                                    z1 z2))
  288.                (else (error error-string z2))))
  289.             ((and operator-ok? (&operator? z1))
  290.              (operation*operator*object->operator operation z1 z2))
  291.             (else
  292.              (error error-string z1))))))
  293.     operation))
  294.  
  295. (define &+
  296.   (make-binary-componentwise-operation r-+ r-+ #t
  297.      &make-rectangular
  298.      "+: bad number"))
  299.  
  300. (define &-
  301.   (make-binary-componentwise-operation r-- r-- #t
  302.      &make-rectangular
  303.      "-: bad number"))
  304.  
  305. (define &=
  306.   (make-binary-componentwise-operation r-= r-= #f
  307.      (lambda (p q) (and p q))
  308.      "=: bad number"))
  309.  
  310. (define (&* z1 z2)
  311.   (cond ((real? z1)
  312.      (cond ((real? z2) (r-* z1 z2))
  313.            ((&complex? z2)
  314.         (&make-rectangular (r-* z1 (&real-part z2))
  315.                    (r-* z1 (&imag-part z2))))
  316.            ((&operator? z2)
  317.         (operation*constant*operator->operator &* z1 z2))
  318.            (else (error "*: bad number" z2))))
  319.     ((&complex? z1)
  320.      (cond ((real? z2)
  321.         (&make-rectangular (r-* (&real-part z1) z2)
  322.                    (r-* (&imag-part z1) z2)))
  323.            ((&complex? z2)
  324.         (&make-rectangular
  325.          (r-- (r-* (&real-part z1) (&real-part z2))
  326.               (r-* (&imag-part z1) (&imag-part z2)))
  327.          (r-+ (r-* (&real-part z1) (&imag-part z2))
  328.               (r-* (&imag-part z1) (&real-part z2)))))
  329.            ((&operator? z2)
  330.         (operation*constant*operator->operator &* z1 z2))
  331.            (else (error "*: bad number" z2))))
  332.     ((&operator? z1)
  333.      (operation*operator*object->operator &* z1 z2))
  334.     (else
  335.      (error "*: bad number" z1))))
  336.  
  337. (define (square z) (&* z z))
  338.  
  339. (define (&/ z1 z2)
  340.   (cond ((real? z2)
  341.      (/real z1 z2))
  342.     ((&complex? z2)
  343.      (/real (&* z1 (&conjugate z2)) (&sqr-magnitude z2)))
  344.     ((&operator? z2)
  345.      (operation*object*operator->operator &/ z1 z2))
  346.     (else
  347.      (error "/: bad number" z2))))
  348.  
  349. (define (/real z x)
  350.   (cond ((r-zero? x)
  351.      (if (fluid-bound? operator-value-retry)
  352.          ((fluid operator-value-retry))
  353.          (error "/: division by zero" z x)))
  354.     ((real? z)
  355.      (r-/ z x))
  356.     ((&complex? z)
  357.      (&make-rectangular (r-/ (&real-part z) x)
  358.                             (r-/ (&imag-part z) x)))
  359.     ((&operator? z)
  360.      (operation*operator*constant->operator /real z x))
  361.     (else
  362.      (error "/: bad number" z))))
  363.  
  364. ;;;; User selector, constructors, and type predicates
  365.  
  366.  
  367. (define complex?
  368.   (lambda (z)
  369.     (or (real? z) (&complex? z))))
  370.  
  371. (define (make-rectangular x y)
  372.   (if (and (real? x) (real? y))
  373.       (&make-rectangular x y)
  374.       (error "Make-rectangular: bad arguments" x y)))
  375.  
  376. (define (make-polar r theta)
  377.   (if (and (real? r) (real? theta))
  378.       (&make-polar r theta)
  379.       (error "Make-polar: bad arguments" r theta)))
  380.  
  381. (define (real-part z)
  382.   (cond ((real? z) z)
  383.     ((&complex? z) (&real-part z))
  384.     ((&operator? z) (operation*operator->operator real-part z))
  385.     (else (error "Real-part: bad number" z))))
  386.  
  387. (define (imag-part z)
  388.   (cond ((real? z) 0)
  389.     ((&complex? z) (&imag-part z))
  390.     ((&operator? z) (operation*operator->operator imag-part z))
  391.     (else (error "Imag-part: bad number" z))))
  392.  
  393. (define (magnitude z)
  394.   (cond ((real? z) (r-abs z))
  395.     ((&complex? z) (&magnitude z))
  396.     ((&operator? z) (operation*operator->operator magnitude z))
  397.     (else (error "Magnitude: bad number" z))))
  398.  
  399. (define-integrable abs magnitude)
  400.  
  401. (define (angle z)
  402.   (cond ((real? z) (if (negative? z) pi 0))
  403.     ((&complex? z) (&angle z (&magnitude z)))
  404.     ((&operator? z) (operation*operator->operator angle z))
  405.     (else (error "Angle: bad number" z))))
  406.  
  407. (define (conjugate z)
  408.   (cond ((real? z) z)
  409.     ((&complex? z) (&conjugate z))
  410.     ((&operator? z) (operation*operator->operator conjugate z))
  411.     (else (error "Conjugate: bad number" z))))
  412.  
  413. ;;;; Unary arithmetic and predicates
  414.  
  415. (define zero?
  416.   (lambda (z)
  417.     (cond ((real? z)
  418.        (r-zero? z))
  419.           ((&complex? z)
  420.        (and (r-zero? (&real-part z))
  421.             (r-zero? (&imag-part z))))
  422.           (else
  423.        (error "Zero?: Bad complex number" z)))))
  424.  
  425. (define-integrable make-unary-componentwise-operation
  426.   (lambda (r-op error-string)
  427.     (letrec ((operation
  428.           (lambda (z)
  429.         (cond ((real? z)
  430.                (r-op z))
  431.               ((&complex? z)
  432.                (&make-rectangular (r-op (&real-part z))
  433.                       (&imag-part z)))
  434.               ((&operator? z)
  435.                (operation*operator->operator operation z))
  436.               (else
  437.                (error error-string z))))))
  438.       operation)))
  439.  
  440. (define 1+
  441.   (make-unary-componentwise-operation r-1+ "1+: bad number"))
  442.  
  443. (define -1+
  444.   (make-unary-componentwise-operation r--1+ "-1+: bad number"))
  445.  
  446. ;;;; N-ary predicates
  447.  
  448. ;;; Predicates are extensions of common Lisp because they can be
  449. ;;; given no arguments, in which case they return `#t'.
  450.  
  451. (define-integrable make-pairwise-test
  452.   (lambda (&pred)
  453.     (lambda args
  454.       (define (loop x y rem)
  455.     (and (&pred x y)
  456.          (or (null? rem)
  457.          (loop y (car rem) (cdr rem)))))
  458.       (or (null? args)
  459.       (null? (cdr args))
  460.       (loop (car args) (cadr args) (cddr args))))))
  461.  
  462. (define = (make-pairwise-test &=))
  463.  
  464. ;;;; N-ary arithmetic
  465.  
  466. (define-integrable accumulation
  467.   (lambda (operation identity)
  468.     (lambda rest
  469.       (define (loop accum rem)
  470.     (if (null? rem)
  471.         accum
  472.         (loop (operation accum (car rem)) (cdr rem))))
  473.       (cond ((null? rest) identity)
  474.         ((null? (cdr rest)) (car rest))
  475.         (else (operation (car rest) (loop (cadr rest) (cddr rest))))))))
  476.  
  477. (define + (accumulation &+ 0))
  478. (define * (accumulation &* 1))
  479.  
  480. (define-integrable inverse-accumulation
  481.   (lambda (operation1 operation2 identity)
  482.     (lambda rest
  483.       (define (loop accum rem)
  484.     (if (null? rem)
  485.         accum
  486.         (loop (operation2 accum (car rem)) (cdr rem))))
  487.       (cond ((null? rest) identity)
  488.         ((null? (cdr rest)) (operation1 identity (car rest)))
  489.         ((null? (cddr rest)) (operation1 (car rest) (cadr rest)))
  490.         (else (operation1 (car rest) (loop (cadr rest) (cddr rest))))))))
  491.  
  492. (define - (inverse-accumulation &- &+ 0))
  493. (define / (inverse-accumulation &/ &* 1))
  494.  
  495. ;;;; Transcendental functions
  496.  
  497. (define (exp z)
  498.   (cond ((real? z) (r-exp z))
  499.     ((&complex? z)
  500.      (&make-polar (r-exp (&real-part z)) (&imag-part z)))
  501.     ((&operator? z)
  502.      (operation*operator->operator exp z))
  503.     (else (error "Exp: bad number" z))))
  504.  
  505. (define (log z)
  506.   (cond ((real? z)
  507.      (if (negative? z)
  508.          (&make-rectangular (error-log (r-- 0 z)) pi)
  509.          (error-log z)))
  510.     ((&complex? z)
  511.      (let ((m (&magnitude z)))
  512.        (&make-rectangular (error-log m) (&angle z m))))
  513.     ((&operator? z)
  514.      (operation*operator->operator log z))
  515.     (else
  516.      (error "log: bad number" z))))
  517.  
  518. (define (error-log x)
  519.   (if (r-zero? x)
  520.       (if (fluid-bound? operator-value-retry)
  521.           ((fluid operator-value-retry))
  522.           (error "log: Log of zero"))
  523.       (r-log x)))
  524.  
  525. (define (expt z1 z2)
  526.   (if (and (real? z1) (real? z2))
  527.       (r-expt z1 z2)
  528.       (exp (&* z2 (log z1)))))
  529.  
  530. (define log10
  531.   (let ((e^base (r-log 10)))
  532.     (lambda (z)
  533.       (/ (log z) e^base))))
  534.  
  535. (define (sqrt z)
  536.   (cond ((real? z)
  537.      (if (negative? z)
  538.          (&make-rectangular 0 (r-sqrt (r-- 0 z)))
  539.          (r-sqrt z)))
  540.     ((&complex? z)
  541.      (let ((m (&magnitude z)))
  542.        (&make-polar (r-sqrt m)
  543.             (r-/ (&angle z m) 2))))
  544.     ((&operator? z)
  545.      (operation*operator->operator sqrt z))
  546.     (else
  547.      (error "sqrt: bad number" z))))
  548.  
  549. (define sin
  550.   (let ((2i (&* 2 i)))
  551.     (lambda (z)
  552.       (cond ((real? z) (r-sin z))
  553.         ((&complex? z)
  554.          (let ((iz (&* i z)))
  555.            (&/ (&- (exp iz) (exp (&- 0 iz)))
  556.            2i)))
  557.         ((&operator? z)
  558.          (operation*operator->operator sin z))
  559.         (else
  560.          (error "sin: bad number" z))))))
  561.  
  562. (define (cos z)
  563.   (cond ((real? z) (r-cos z))
  564.     ((&complex? z)
  565.      (let ((iz (&* i z)))
  566.        (&/ (&+ (exp iz) (exp (&- 0 iz)))
  567.         2)))
  568.     ((&operator? z)
  569.      (operation*operator->operator cos z))
  570.     (else
  571.      (error "cos: bad number" z))))
  572.  
  573. (define (tan z)
  574.   (cond ((real? z)
  575.      (let ((den (r-cos z)))
  576.        (if (r-zero? den)
  577.                (if (fluid-bound? operator-value-retry)
  578.                    ((fluid operator-value-retry))
  579.                (error "tan: Overflow" z))
  580.            (r-/ (r-sin z) den))))
  581.     ((&complex? z)
  582.      (let* ((iz (&* i z))
  583.         (den (&+ (exp iz) (exp (&- 0 iz)))))
  584.        (if (zero? den)
  585.            (if (fluid-bound? operator-value-retry)
  586.                    ((fluid operator-value-retry))
  587.                (error "tan: Overflow" z))
  588.             (&* -i (&/ (&- (exp iz) (exp (&- 0 iz))) den)))))
  589.     ((&operator? z)
  590.      (operation*operator->operator tan z))
  591.     (else
  592.      (error "tan: bad number" z))))
  593.  
  594. (define atan 
  595.   (let ((2i (&* 2 i)))
  596.     (lambda (y . optionals)
  597.       (let ((x (if (not (null? optionals)) (car optionals) 1)))
  598.         (if (and (real? y) (real? x))
  599.             (r-atan y x)
  600.             (let ((iy (&* i y)))
  601.               (&/ (&- (log (&+ x iy)) (log (&- x iy))) 2i)))))))
  602.  
  603. (define (acos z)
  604.   (if (and (real? z) (<= z 1) (>= z -1))
  605.       (r-acos z)
  606.       (&* -i (log (&- z (&* -i (sqrt (&- 1 (&* z z)))))))))
  607.  
  608. (define (asin z)
  609.   (if (and (real? z) (<= z 1) (>= z -1))
  610.       (r-asin z)
  611.       (&* -i (log (&+ (&* i z) (sqrt (&- 1 (&* z z))))))))
  612.  
  613. (define-integrable number? complex?)
  614.  
  615. (define operator? 
  616.   (lambda (z)
  617.     (or (number? z) (&operator? z))))
  618.  
  619.  
  620. ;;; Because there is no rationals here.
  621.  
  622. (define make-rational /)
  623.  
  624.  
  625. ;;; 6.003 useful operators
  626.  
  627. (define s identity-operator)
  628. (define t identity-operator)
  629.  
  630. ;;; For electrical engineers
  631.  
  632. (define j  i)
  633. (define -j -i)
  634.